

function  [mu_t0, Sigma_t0] = gp_update(ya, ya_hist, eta_hist, sigma_etasq_hist, mu0,sigma_c,psi, cov_fun)

%  inputs
%  parameters: r, ybar, rho_y, W_cc, kappa

 
%  prior on chat(y) and Sigma(y,y')
if isempty(ya_hist)

    Sigma_t0 = sigma_c^2*ones(size(ya)); %Var(cstar(y) | eta_hist )
    mu_t0 = mu0(ya);
        
else
    
    if cov_fun == 1
        temp1 = bsxfun(@minus, ya, ya_hist');
        %temp2 = bsxfun(@minus, a_t, a_hist');
        Sigma12  = sigma_c^2*exp(-psi*(temp1.^2));
        
        temp1 = bsxfun(@minus, ya_hist, ya_hist');
        %temp2 = bsxfun(@minus, a_hist, a_hist');
        Sigma_y1yN = sigma_c^2*exp(-psi*(temp1.^2));
        Sigma22 = Sigma_y1yN + diag(ones(length(ya_hist),1).*(sigma_etasq_hist));
        
        L = chol(Sigma22, 'lower');
        v = L\Sigma12';
        %Sigma_t0 = diag(sigma_c^2 - v'*v);  %Var(cstar(y) | eta_hist )
        
        sum_vsq = sum(v.^2,1);
        Sigma_t0 = (sigma_c^2 - sum_vsq)';
        
        
        mu_t0 = (L')\(L\(eta_hist - mu0(ya_hist )));
        mu_t0 = mu0(ya) + Sigma12*mu_t0;
    else
        temp1 = bsxfun(@minus, ya, ya_hist');
        %temp2 = bsxfun(@minus, a_t, a_hist');
        Sigma12  = sigma_c^2*exp(-psi*(abs(temp1)));
        
        temp1 = bsxfun(@minus, ya_hist, ya_hist');
        %temp2 = bsxfun(@minus, a_hist, a_hist');
        Sigma_y1yN = sigma_c^2*exp(-psi*(abs(temp1)));
        Sigma22 = Sigma_y1yN + diag(ones(length(ya_hist),1).*(sigma_etasq_hist));
            
        L = chol(Sigma22, 'lower');
        v = L\Sigma12';
        Sigma_t0 = diag(sigma_c^2 - v'*v);  %Var(cstar(y) | eta_hist )
            
        mu_t0 = (L')\(L\(eta_hist - mu0(y_hist + a_hist)));
        mu_t0 = mu0(ya) + Sigma12*mu_t0;
    end
    
end
    


end